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Radioactive isotopes produced through cosmic muon spallation are a background for rare-event detection in 
v detectors, double-/3-decay experiments, and dark-matter searches. Understanding the nature of cosmogenic 
backgrounds is particularly important for future experiments aiming to determine the pep and CNO solar neu- 
trino fluxes, for which the background is dominated by the spallation production of C. Data from the Kamioka 
liquid scintillator antineutrino detector (KamLAND) provides valuable information for better understanding 
these backgrounds, especially in liquid scintillators, and for checking estimates from current simulations based 
upon MUSIC, FLUKA, and GEANT4. Using the time correlation between detected muons and neutron captures, 
the neutron production yield in the KamLAND liquid scintillator is measured to be Y n = (2.8 ± 0.3) x 10 -4 
^i _1 g _1 cm 2 . For other isotopes, the production yield is determined from the observed time correlation related 
to known isotope lifetimes. We find some yields are inconsistent with extrapolations based on an accelerator 
muon beam experiment. 
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I. INTRODUCTION 

Cosmic-ray muons and their spallation products are poten- 
tial sources of background for neutrino detectors, double-/?- 
decay experiments, and dark-matter searches, even when the 
detectors are deployed underground. Characterizing cosmic- 
ray-muon-induced backgrounds, particularly the secondary 
neutrons and radioactive isotopes produced by muon-initiated 
spallation processes, is essential for interpreting these experi- 
ments. 

Liquid-scintillator detectors such as the KamLAND 
(Kamioka liquid-scintillator antineutrino detector), Borex- 
ino BjJ— 13|] , CANDLES IV (calcium fluoride for studies of neu- 
trino and dark matters by low-energy spectrometer) [4, 5], 
SNO+ (Sudbury Neutrino Observatory) H 0], LENS (low 
energy neutrino spectroscopy) [8], and LENA (low-energy 
neutrino astronomy) [9] are designed to detect low-energy 
phenomena. In an organic liquid scintillator (LS), energetic 
muons and subsequent showers interact mostly with 12 C, the 
most abundant nucleus heavier than 1 H in the LS, generat- 
ing neutrons and isotopes by electromagnetic or hadronic pro- 
cesses. The muon-initiated spallation of carbon targets is a 
matter of primary interest. 

Isotope production by muon-initiated spallation has been 
studied by an earlier experiment [ 10] using the CERN Super 
Proton Synchrotron (SPS) muon beam. The energy depen- 
dence was studied with 100 and 190 GeV incident muons. 
The production yield at other energies is estimated from these 
data by extrapolation, assuming a power-law dependence on 
the muon energy. Direct measurements of the production 
yield by underground detectors such as the Mont Blanc liq- 
uid scintillation detector (LSD) 111 111 , the INFN large-volume 
detector (LVD) [12], and Borexino [13] were compared with 
calculations exp loiting simulations based on MUSIC Q, 
FLUKA 111 Gil, and GEANT4 (nl [Hi. Particular attention 
was paid to neutron production since isotope production mea- 
surements are difficult with the small scintillator masses used 
in these detectors. KamLAND, owing to its larger mass — ~1 
kiloton of LS — does not suffer from this difficulty and is well 
placed to study a variety of isotopes of interest. 

This paper presents the neutron and isotope production 
rates in KamLAND from muon-initiated spallation based 
upon data collected from 5 March 2002 to 12 May 2007. The 
results are compared with simulations and other experiments. 
These comparisons provide important information for validat- 
ing Monte Carlo simulations. 



II. DETECTOR DESCRIPTION AND PERFORMANCE 

KamLAND is located under the peak of Ikenoyama (Ike 
Mountain, 36.42°N, 137.31°E), and the vertical rock over- 
burden is approximately 2700 meters water equivalent (mwe). 
A schematic diagram of KamLAND is shown in Fig. Q] 
KamLAND consists of an active detector region of approxi- 
mately 1 kiloton (kton) of ultrapure LS contained in a 13-m- 
diameter spherical balloon made of 135-/im-thick transparent 
nylon-EVOH (ethylene vinyl alcohol copolymer) composite 




FIG. 1: Schematic diagram of the KamLAND detector at the 
Kamioka underground laboratory. 



TABLE I: Elemental composition of the KamLAND liquid scintil- 
lator. The hydrogen-to-carbon ratio was verified by elemental anal- 
ysis (±2%). The liquid scintillator contains traces of nitrogen and 
oxygen from its fluor, PPO (C15H11NO), and dissolved gases. Mea- 
surements in similar liquid hydrocarbons 1 19, 20] indicate that the 
amount of dissolved nitrogen saturates at about 200-500 parts per 
million (ppm). The range of nitrogen composition values given in 
this table reflects the extreme cases of zero dissolved nitrogen gas 
to full saturation. The dissolved oxygen content of liquid scintillator 
taken from the center of KamLAND was measured to be less than 



3 ppm, which is insignificant compared to the oxygen contribution 
from PPO. 







Number of targets 


Element 


Stoichiometry 


(per kiloton) 


Hydrogen 


1.97 


8.47 x 10 31 


Carbon 


= 1 


4.30 x 10 31 


Nitrogen 


1 x 10~ 4 to 6 x 10" 4 


5 x 10 27 to 3 x 10 28 


Oxygen 


1 x 10~ 4 


5 x 10 27 



film and supported by a network of Kevlar ropes. In addi- 
tion to providing containment for the LS, the balloon protects 
the LS against the diffusion of ambient radon from the sur- 
rounding components. The total volume of LS in the bal- 
loon is 1171 ± 25 m 3 as determined with flow meters dur- 
ing the initial filling. The LS consists of 80% dodecane and 
20% pseudocumene (l,2,4-trimethylbenzene)by volume, and 
1.36 ±0.03 g/liter of the fluor PPO (2,5-diphenyloxazole). 
The density of the LS is 0.780 g/cm 3 at 11.5°C. The calcu- 
lated elemental composition of the LS is given in Table|T] 

A buffer comprising 57% isoparaffin and 43% dodecane 
oils by volume fills the region between the balloon and the sur- 
rounding 18-m-diameter spherical stainless-steel outer vessel 
to shield the LS from external radiation. The specific gravity 
of the buffer oil (BO) is adjusted to be 0.04% lower than that 
of the LS. An array of photomultiplier tubes (PMTs) — 1325 
specially developed fast PMTs masked to 17-in. -diameter and 
554 older 20-in. -diameter PMTs reused from the Kamiokande 



experiment [21] — are mounted on the inner surface of the 
outer containment vessel, providing 34% photocathode cov- 
erage. During the period from 5 March 2002 to 27 Febru- 
ary 2003 the photocathode coverage was only 22%, since the 
20-in. PMTs were not operated. A 3-mm-thick acrylic bar- 
rier at 16.6 m in diameter helps prevent radon emanating from 
the PMT glass from entering the BO. The inner detector (ID), 
consisting of the LS and BO regions, is surrounded by a 3.2- 
kton water-Cherenkov detector instrumented with 225 20-in. 
PMTs. This outer detector (OD) absorbs 7 rays and neutrons 
from the surrounding rock and enables tagging of cosmic-ray 
muons. 

The KamLAND front-end electronics (FEE) system 
is based on the analog transient waveform digitizer 
(ATWD) JH which captures PMT signals in 128 10-bit digi- 
tal samples at intervals of 1.5 ns. Each ATWD captures three 
gain levels of a PMT signal to obtain a dynamic range from 
1 photoelectron (p.e.) to 1000 p.e. Each ATWD takes 27 /xs 
to read out, so two are attached to each PMT channel to re- 
duce dead time. The FEE system contains discriminators set 
at 0.15 p.e. (^0.3 mV) threshold which send a 125-ns-long 
logic signal to the trigger electronics. The trigger electronics 
counts the number of ID and OD PMTs above the discrimi- 
nator threshold with a sampling rate of 40 MHz and initiates 
readout when the number of 17-in. ID PMTs above the dis- 
criminator threshold (A17) exceeds the number correspond- 
ing to ^0.8 MeV deposited energy. The trigger system also 
issues independent readout commands when the number of 
OD PMTs above threshold exceeds a preset number. 

The energy can be estimated from the A max parameter, de- 
fined as the maximum value of A17 in a 200-ns period follow- 
ing the trigger command. However, the offline analysis takes 
full advantage of the information stored in the digitized PMT 
signals by identifying individual PMT pulses in the waveform 
information that is read out. The time and integrated area 
(called "charge") are computed from the individual pulses. 
For each PMT, the average charge corresponding to a single 
p.e. is determined from single-pulse waveforms observed in 
low-occupancy events. The ID PMT timing is calibrated with 
light pulses from a dye laser (^1.2-ns pulse width), injected 
at the center of the detector through an optical fiber. The ver- 
tices of spatially localized low-energy (<30 MeV) events are 
estimated by comparing calculated time-of-flights of optical 
photons from the hypothetical vertex to the measured arrival 
times at the PMTs in KamLAND. 

The reconstructed energies of events are calibrated with 7 
sources: 203 Hg, 68 Ge, 65 Zn, and 60 Co; and with n+7 sources: 
241 Am+ 9 Be and 210 Po+ 13 C JH. These are deployed at var- 
ious positions along the vertical axis of the detector and oc- 
casionally off the vertical axis within 5.5 m from the detector 
center [24]. Such calibrations cover energies between 0.28 
and 6.1 MeV. The energy calibration is aided with studies 
of background contaminants 40 K and 208 T1, 212 Bi- 212 Po and 
214 Bi- 214 Po sequential decays, 12 B and 12 N spallation prod- 
ucts, and 7's from thermal neutron captures on 4 H and 12 C. 

The visible energy E v i s of an event is computed from the 
measured light yield. Specifically, _E V j S is the number of de- 
tected p.e. after corrections for PMT variation, dark noise, 
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FIG. 2: Ikenoyama topological profile J27h . The black point near the 
center is the location of KamLAND. 



solid angle, shadowing by suspension ropes, optical trans- 
parencies, and scattering properties in the LS. The relation- 
ship between _E V ; S and the deposited energy £dep of 7's, e ± 's, 
protons, and a's is nonlinear and modeled as a combination 
of Birks-quenched scintillation 1I25I l26Tl and Cherenkov radi- 
ation. The scale is adjusted so that i5 v j s is equal to i?dep for 
the 2.225-MeV 7 ray from neutron capture on 1 H. The ob- 
served energy resolution is ^7.4%/ -\/i? V j S (MeV) for the pe- 
riod without the 20-in. PMTs, and ~6.5%/y/E vis (MeV) for 
the rest of the data. 

The calibration sources are also used to determine system- 
atic deviations in position reconstruction by comparison with 
the source's known position. This comparison gives an aver- 
age position reconstruction uncertainty of less than 3 cm for 
events with energies in the range 0.28 to 6. 1 MeV. 



III. COSMIC-RAY MUONS 

A digital map [27] of the topological profile of Ikenoyama 
is shown in Fig. [2] The vertical overburden at KamLAND is 
approximately 1000 m of rock and the minimum overburden 
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FIG. 3: Correlation between the light yield in the inner detector and 
the shortest distance from the muon track to the KamLAND detector 
center (impact parameter). The vertical dashed red line represents 
the boundary between the LS and BO at 650 cm. Events located to 
the left of this boundary and above the horizontal blue dashed line 
(Cm = 4 x 10 4 p.e.) are designated as LS muons. 
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corresponding to a nearby valley is approximately 900 m. 

Cosmic-ray muons are identified either by the large amount 
of scintillation and Cherenkov light detected by the ID PMTs, 
or by the Cherenkov light detected by the OD PMTs. Muons 
crossing the ID (ID muons) are selected by requiring that one 
of the following conditions was satisfied: 

1. Ad > 10000 p.e. (~30MeV), 

2. Ad > 500 p.e. and N OD > 5, 

where Ad is the total light yield measured by the ID 17-in. 
PMTs, and A*od is the number of OD PMTs with signals 
above threshold. Approximately 93% of the ID muons sat- 
isfy the first selection criterion. The ID muon track is re- 
constructed from arrival times of the first-arriving Cherenkov 
or scintillation photons at the PMTs. Since for relativistic 
muons the wavefront of the scintillation light proceeds at the 
Cherenkov angle, and since muons generate enough light to 
generate photoelectrons in every PMT, by restricting the fit 
to the first-arriving photons both Cherenkov and scintillation 
photons can be treated identically. The observed muon track is 
then established by minimizing time-of-flight deviations from 
hypothetical muon tracks. The fit converges for 97% of all 
ID muon events. The majority of the events that are not re- 
constructed are believed to be multiple muons or muons ac- 
companied by large electromagnetic or hadronic showers for 
which the tracking model is not valid. 

Muons passing only through the BO produce mostly 
Cherenkov light, whereas muons passing through the LS gen- 
erate both Cherenkov and scintillation light. Figure [3] shows 
the correlation between the light yield (Ad) and the short- 
est distance between the reconstructed muon track and the 
center of KamLAND (impact parameter). The boundary at 
650 cm between the BO and LS regions is evident. The corre- 
lations between Ad and the reconstructed muon track length 
in the BO and LS regions (Lbo and Lls, respectively) are 



FIG. 4: Correlation between the total light yield measured by the 
ID 17-in. PMTs (£id) and the muon track length for (a) muons 
that pass through both the LS and BO regions and (b) muons that 
pass through only the BO region. The red lines show the fitted 
light yield per unit length of (dCs/dX) — 629 ± 47 p.e./cm and 
(dC 6 /dX) = 31 ±2 p.e./cm. 

plotted in Figs.Ua) and@|b). A linear trend, corresponding to 
minimum ionizing muons, is apparent in both distributions. 
The slope of each line is the light yield per unit length in 
the respective material. In BO, where the light is predomi- 
nantly Cherenkov, the light yield per unit length is found to 
be (dCg/dX) — 31 ± 2 p.e./cm; the fit was restricted to path 
lengths above 700 cm, since fits at shorter path lengths are 
complicated by the presence of PMTs which may obstruct 
some of the emitted light. In the LS we obtain 

= ClU L ™ (dCc/dX) = 629 ± 47 p.e./cm, (1) 

where dCs/dX includes the Cherenkov light created in the 
LS. The muons in Fig. [4] generating light yields above the 
baseline linear trend are likely to involve secondary particles. 
We define an excess light yield parameter AC, 

< 2 > 

for the purpose of describing showering muons associated 
with secondary particles. 

The LS muon rate is estimated by selecting muons with 
Ad > 4 x 10 4 p.e. and impact parameter < 650 cm. The 
light yield cut has a negligible inefficiency for LS muons, 
while the impact parameter cut eliminates LS muons that 
reconstruct outside the balloon because of the resolution of 
the fitter algorithm. This provides the lower limit on the 
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TABLE II: Calculated J M and E^ values for various rock types 
using the simulation method from Ref. 13311 . The range in val- 
ues for J M and corresponds to varying the specific gravity of 
the rock from 2.65 to 2.75 g/cm 3 . The generic skarn is defined 
to be 70% granite and 30% calcite (by weight) with the following 
properties: <Z>= 10.22, <A>= 20.55 a.u., and radiation length 
X = 25.411 g/cm 2 . 



Ikenoyama rock model 


J M (m- 2 h" 1 ) 


E» (GeV) 


Inishi rock 


5.66-6.71 


262-268 


Standard rock 


4.95-5.83 


256-262 


Generic skarn 


4.90-5.82 


254-260 


This measurement 


5.37 ±0.41 





FIG. 5: Distribution of muon track lengths in the KamLAND LS. The 
dashed blue line shows the expected linear distribution for a sphere of 
nominal volume, and the average track length in this case is 867 cm. 
The average muon track length measured in the detector is 878 cm, 
in agreement with the value 874 ± 13 cm calculated using the proper 
LS shape and muon angular distributions. 



LS muon rate. An upper limit is established by remov- 
ing the impact parameter cut and increasing the £id cut to 
> 10 5 p.e. to eliminate muons that pass through the BO with- 
out transversing the LS. This cut again has a small ineffi- 
ciency for LS muons but does not eliminate muons which 
shower in the BO. Averaging the two measurements yields 
= 0.198 ± 0.014 Hz, where the error estimate is domi- 
nated by systematic uncertainties due to temporal variations, 
the difference between the two measurements, and the pres- 
ence of muon bundles, in which multiple simultaneous muons 
may be tagged as single-muon events. Using the parametriza- 
tion of Ref. [28], the effect of muon bundles was estimated to 
be <5%. The muon rate corresponds to an integrated muon 
intensity of = 5.37 ± 0.41 m~ 2 h _1 , where the error in- 
cludes uncertainties in the shape of the LS volume but is dom- 
inated by the uncertainty in the muon rate. Using a cut of 
AC > 10 6 p.e., which is equivalent to a ^3 GeV threshold, 
the rate of showering muons in the LS is ^0.03 Hz. It is pos- 
sible that some atmospheric neutrino interactions leak into the 
LS muon sample. However, an estimate with the neutrino flux 
from Ref. B29I1 and cross sections from Ref. [30] gives less 
than 4 x 10 -5 Hz from atmospheric neutrinos. The muon 
track length distribution is shown in Fig. [5] The measured 
average track length is = 878 cm, in agreement with the 
calculated value of L M = 874 ± 13 cm where the nonspher- 
ical corrections to the balloon shape and the muon angular 
distributions were taken into account. 

The production yield of radioactive isotopes from muon- 
initiated spallation is energy dependent. KamLAND does not 
measure the muon energy, so it is estimated from simulation. 
There are previous estimates of the mean energy ranging 
from 198 GeV (H to 285 GeV Gil HI at the Kamiokande, 
KamLAND, and SuperKamiokande detectors. The authors of 
Ref. A3 311 made a detailed determination of J M and using 
the muon simulation code (MUSIC) ||14tl to transport muons, 
generated according to the modified Gaisser parameteriza- 



tion 113411 sea-level muon flux distribution, through a digital 
profile of Ikenoyama. Although this calculation reproduces 
the zenith and azimuthal angular distributions observed by 
KamLAND, as shown in Fig. 10 of Ref. [33], it overestimates 
J M by 14% relative to this work. The calculation assumed a 
homogeneous rock entirely of the Inishi type (see Table II of 
Ref. 13311 for the chemical composition), but other rock types, 
such as granite, limestone, and several types of metamorphic 
rock, are common in Ikenoyama in unknown quantities 13511 . 
Table [II] shows the result of calculations of J M and E^ us- 



ing the same simulation method of Ref. 113311 but for stan- 
dard rock lf36L l37ll and generic skarn, a generic mixture of 
rock types found in a skarn-type mine like that of Ikenoyama. 
Here, generic skarn is defined to be 70% granite and 30% cal- 
cite by weight. The J M values range from 4.90 m~ 2 h _1 , for 
2.75 g/cm 3 specific gravity generic skarn, to 6.71 m -2 h _1 , 
for 2.65 g/cm 3 specific gravity Inishi rock, while varies 



in excellent 



from 254 GeV to 268 GeV. The value of J M for 2.70 g/crrC 
specific gravity standard rock is 5.38 m~ 2 h -1 , 
agreement with our measured value. The value of Eu for this 
rock is 259 GeV. We take E^ = 260 ± 8 GeV, where the un- 
certainty is chosen to cover the full range for the various rock 
types. 



IV. SPALLATION NEUTRON YIELD 

Most of the neutrons produced in the KamLAND LS 
are captured on hydrogen or carbon atoms. The capture 
cross section varies inversely with respect to velocity, and 
the mean neutron capture time r n is constant with respect 
to energy. The capture time (t) distribution is exponen- 
tial, P(t) oc e~*/ r ". A calculation using the elemental com- 
position of the KamLAND LS given in Table U and the 
thermal neutron capture cross sections from Ref. [38] gives 
t„ = 206 [is. This calculation indicates that 99.5% of the neu- 
trons are captured on 1 H, while the remainder are captured 
mostly on 12 C. The probability for capture on the other iso- 
topes in the LS, such as 13 C, is 2 x 10 -4 or less. 

Neutrons produced by muon-initiated spallation in the LS 
can be identified by the characteristic capture 7 rays. Figure|6] 
shows the _E V i S distributions in signal (150 < AT < 1000/is) 
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FIG. 6: Visible energy spectra of the events following muons in the 
signal time window 150 < AT < 1000/is (crosses) compared with 
background in the time window 4150 < AT < 5000/is (dashes). 
The low-energy tails in the n}\{ and n 12 C peaks are caused by 
the electronics effect discussed in the text. In the time window 
1300 < AT < 2150/^s (solid) there is a small n 1 !! peak without the 
low-energy tail. 

and background (4150 < AT < 5000/xs) coincidence win- 
dows following some muons, where AT = (i — t^) is the 
time elapsed since the muon's passage. The E V - 1S distri- 
bution clearly shows peaks from neutron captures on : H 
(2.225 MeV) and 12 C (4.9 MeV), which are not evident in 
the background. Figure [7] shows the AT distribution for 
events within the LS volume and with 1.8 < E v \ s < 2.6 MeV, 
which includes the single 2.225-MeV 7 ray emitted by neu- 
tron capture on 1 H. For AT < 1000/zs, there is a clear devia- 
tion from the exponential distribution due to the overload that 
large muon signals produce on individual electronics channels 
and to the dead time in the system arising from the very high 
event multiplicity following the muon. Both effects intervene 
in events that are quite different from those that KamLAND 
was designed to record. 

The number of neutrons produced by muon-initiated spalla- 
tion in the LS is established by a binned maximum likelihood 
fit 113911 to the data in Fig. [TJ using the function 

r (t) = ^ e -(*-**)/T» +rBj (3 ) 

where N n is the total number of neutron captures associated 
with the selected events, r„ is the mean neutron capture time, 
and tb is the background rate, which is assumed to be ap- 
proximately constant in the region of interest (AT < 2500/is) 
because of the low muon rate (~ 0.2 Hz). To avoid the 
electronics-induced distortions, the fit is restricted to the re- 
gion AT > 1300/zs. The parameters N n and tb are free, but 
the mean capture time is constrained to r„ = 207.5 ± 2.8/xs 
with a Gaussian penalty function. This mean capture time 
is determined from two independent measurements of r„: a 
241 Am+ 9 Be calibration source, and an analysis of a sample 
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AT (lis) 

FIG. 7: Time difference between a muon and the neutron capture 
event, with the cuts 1.8 < E v \s < 2.6 MeV. The red line shows 
the fit restricted to the region AT > 1300/is. The fit results in 
N n = (4.2 ± 0.3) x 10 6 with x 2 / d -°- f ■= 98/118. 



TABLE III: Summary of the dominant contributions to the neutron 
detection efficiency. 



Effect 


Value 


Neutron-eliminating reactions, e.g., (n,p) 


(96.3 ±3.7)% 


Neutron captures on 1 H 


(99.5 ±0.1)% 


LS-BO boundary 


(93.3 ± 2.0)% 


Electronics dead time effects 


> 98% 


Combined efficiency 


(89.4 ± 3.8)% 



of neutrons generated by clean muons. These clean muons 
are identified by a multiplicity parameter t/t defined to be 
the number of trigger commands that follow the muon within 
a 10-ms period. Fits to the AT distribution of the subset 
of neutrons selected with various limits on i]t demonstrate 
that muons with i]t < 10 give unbiased fit residuals. The 
value of r„ is 207.5 ± 0.3/xs from the clean muon sample 
and 205.2 ± 0.5/zs from the 241 Am+ 9 Be source data. The ob- 
served 2.3 /us discrepancy between these values is not com- 
pletely understood but is suspected to be caused by neu- 
trons from the 241 Am+ 9 Be source that are captured on the 
stainless-steel source capsule. In this analysis, we use the 
value r n — 207.5/is from the nonshowering muon sample 
with an uncertainty of ±2.8/is covering both measurements. 
The fit shown in Fig.|7]in the region AT > 1300/is results in 
N n = (4.2 ± 0.3) x 10 6 and x 2 /d.o.f.= 98/118. 

The actual number of neutrons Af n produced by muon- 
initiated spallation is related to the fit result N n by an effi- 
ciency e n , 

K = — , (4) 

En 

which accounts for other neutron-eliminating nuclear reac- 
tions such as 12 C(n,p) and a net neutron loss at the LS-BO 
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boundary. This efficiency is calculated using the MUSlC-based 
muon simulation described in Sec.|III]and the GEANT4-based 
Monte Carlo of KamLAND described in Sec. lVI Al The muon 
simulation generates muons with a three-momentum distri- 
bution appropriate to the KamLAND site inside Ikenoyama. 
Muon transport is then modeled through KamLAND. In this 
simulation neutrons are created and destroyed; neutrons that 
survive to thermalization are tracked until they are captured. 
Energy scale nonlinearities and the finite 7 ray resolution are 
included. Two competing effects of the LS-BO boundary are 
taken into account. The first is for neutrons produced by 
muons in the LS that leak out, leading to an under count- 
ing of neutrons. The second is for neutrons produced outside 
of the LS that leak in, primarily from the BO, leading to an 
over count. The leak-out fraction [(11.7 ± 1.9)%] is larger 
than the leak-in fraction [(5.4 ±0.4)%] resulting in a net 
(6.7 ± 2.0)% correction. The component of the efficiency re- 
lated to the electronics dead time for high-multiplicity events 
is measured by comparing the number of recorded wave- 
forms with the number of trigger commands. The correc- 
tion is measured to be less than 2% for AT > 1300/is. In 
counting neutrons, we use the definition by which the (n, 2n) 
reaction generates one new neutron and the (n, n') reaction 
generates no new neutrons. This method gives an efficiency 
e„ = (89.4 ± 3.8)% (broken down in Table|III). 

Using Eq. © we find Af n = (4.7 ± 0.4) x 10 6 , which is 
then used to extract the neutron production yield: 
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FIG. 8: Background-subtracted T v i s spectrum above 4 MeV. 
Signal and background events are taken from 2 < AT < 60 ms 
and 502 < AT < 560 ms, respectively. The production rate of 
12 B (r = 29.1 ms, Q = 13.4 MeV) is estimated to be 54.8 ± 1.5 
kton _1 day~ , from fitting Eq. (JSJ to the AT distribution shown in 
the inset. This fit has a x 2 /d.o.f.= 509/495. The production rate 
12 N (r = 15.9 ms, Q = 17.3 MeV) is estimated to be 2.2 ± 0.5 
kton _1 day~ , from a similar fit to the AT distribution for events 
with 14 < Tvis < 20 MeV, where the higher E v i a threshold is im- 
posed to exclude 12 B. 



Y n = 



Af n 



(5) 



where R 
LS 



muons, 



= 0.198 ±0.014 Hz is 
T L = 1.24 x 10 s s is 



the measured rate of 
the detector live time, 



p = 0.780 ± 0.001 g/cm d is the density of the LS, and 
L M = 874 ± 13 cm is the calculated mean muon track 
length. The resulting neutron production yield is 
Y n = (2.8 ± 0.3) x 10~ 4 fi^g^cm 2 , with (64 ± 5)% of 
the neutrons produced by events classified as showering 
muons. 



V. SPALLATION ISOTOPE YIELD 



rate of uncorrected background events that are in accidental 
coincidence with muons. 

In terms of Afi, the spallation production yield for isotope i 
is equal to 



Y 



(1) 



where T^, p, and are defined as in Eq. (0. The spal- 
lation production rate for isotope i is equal to 



R, 



pV T T L " 



(8) 



where Vt — 1171 ± 25 m 3 is the target volume. 



The method for determining the yields of spallation- 
generated isotopes is similar to the neutron analysis described 
in Sec. [IV] Spallation-generated isotopes are identified by 
their decay time relative to their creation and by their decay 
energy. The decay time, AT = t — t^, is calculated for each 
event relative to all previous muons. Usually, several different 
isotopes decay in a given time window. The number of each 
isotope produced, Afi, is obtained from a binned maximum 
likelihood fit l39ll to the AT distribution, using the function 

K*) = £— e- (t ^ )/r '+r B , (6) 

1 

where Ni is related to Afi by an event selection efficiency 
(Afi = Ni/ei), Ti is the mean lifetime, and tb is a constant 



A. 12 Banrf 12 N 

12 B (r = 29.1ms, Q = 13.4 MeV)j40ll 0~ decay and 
12 N (r = 15.9 ms, Q = 17.3 MeV) CT/^ decay candi- 
date events are selected via cuts on E v i s and AT. The in- 
set in Fig. [8] shows the distribution of AT for events with 
4 < E vis < 20 MeV and 2 < AT < 500 ms. A binned max- 
imum likelihood fit to the AT distribution using Eq. (O, and 
including the long-lived isotopes 8 He, 9 Li, 9 C, and 12 N as 
contaminants, yields A^( 12 B) = (5.94 ±0.10) x 10 4 events 
with x 2 /d.o.f.= 509/495. Longer lived isotopes give roughly 
constant decay rates on this time scale and fit out as a com- 
ponent of tb- A similar fit to the AT distribution for events 
with 14 < T vis < 20 MeV, where the higher T vis threshold is 
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FIG. 9: Background-subtracted T v i s spectrum above 4 MeV, where 
signal and background events are taken from 0.6 < AT < 4.0 s 
and 10.6 < AT < 14.0 s, respectively. The production rates of 8 Li 
(r = 1.21 s,Q = 16.0 MeV) and 8 B (r = 1.11 s, Q = 18.0 MeV) 
are estimated to be 15.6 ± 3.2 and 10.7 ± 2.9 kton _1 day _1 , respec- 
tively, from simultaneously fitting the _B v is spectrum and the AT 
distribution shown in the inset. The fit to the AT distribution has a 
X 2 /d.o.f- 95/91. 

imposed to exclude 12 B, gives N( 12 N) = (2.8 ± 0.3) x 10 2 
events. A comparison with the predicted E V - 1S spectrum is 
shown in Fig. [8] The fJ V i S spectra are predicted from the 
allowed 12 B and 12 N /3 ± -decay spectra, taking into account 
the KamLAND detector response, normalized to the observed 
iV( 12 B) and N( 12 N). The detector response model includes 
the energy nonlinearities and boundary effects described in 
Sec.lII] 

The inefficiency in identifying 12 B and 12 N candidates is 
dominated by the E v [ s cut. The efficiencies calculated by in- 
tegrating the predicted E V [ S spectra over the selection window 
give e( 12 B) = (82.9 ± 0.7)% and e( 12 N) = (9.3 ± 1.6)%, 
where the errors come from the uncertainty in the detector 
response. Using Eq. the resultant isotope production 
yields are calculated to be F( 12 B) = (42.9 ± 3.3) x 10~ 7 
andF( 12 N) = (1.8 ±0.4) x 10~ 7 ^g^cm 2 . The produc- 
tion rates are i?( 12 B) = 54.8 ± 1.5 and R( 12 N) = 2.2 ± 0.5 
kton^ 1 day _1 . 



B. 8 Li and 8 B 

8 Li (r = 1.21 s, Q = 16.0 MeV) El /3" decay and 8 B 
(r = 1.11 s, Q = 18.0 MeV) HH /3 + decay candidate events 
are selected according to E V i S and AT. The inset in 
Fig. [9] shows the distribution of AT for all preceding muons 
for events with 4 < E vis < 20 MeV and 0.6 < AT < 10 s, 
where the lower limit on AT is chosen to exclude 9 C, 
8 He, 9 Li, and other isotopes with shorter lifetimes. Iso- 
topes with longer lifetimes are roughly constant over the se- 
lected range of AT. To avoid a large accidental coinci- 




Distance from Muon Track (cm) 

FIG. 10: Impact parameter distribution for events identified as 12 B 
(blue points) and neutrons (black histogram) produced by muon- 
induced spallation. Only nonshowering muons (A£ < 10 6 p.e.) 
were used to make this plot. The efficiencies for the AL < 3 m 
track cut (represented by the vertical dashed line) are evaluated to be 
(91.6 ± 4.3)% from 12 B events, and (95.9 ± 0.7)% from neutron 
events. 



dence background from uncorrected muons, all showering 
muons and a portion of nonshowering muons whose track 
is within 3 m of a 8 Li or 8 B candidate are selected. Fig- 
ure [9] shows the E V [ S distribution of 8 Li and 8 B candidate 
events with 0.6 < AT < 4.0 s after subtracting background 
estimated from the range 10.6 < AT < 14.0 s. iV( 8 Li) and 
A^( 8 B) are determined from a simultaneous binned maxi- 
mum likelihood fit to the AT distribution and a chi-square 
fit to the _E V i S distribution. The expected E V i S spectra for 
KamLAND are calculated by convolving the /3 spectra from 
Refs. JUIH] with KamLAND's detector response. For the fit 
to the E v i s distribution, the energy scale parameters are con- 
strained to an allowed region determined by a prior fit to 7 ray 
calibration data (described in Sec. [IT]) and the 12 B T V j S distri- 
bution in Fig. [8] 

The efficiency is given by the product of eg for the 
E V [ S selection and 65 for the muon- 8 Li or muon- 8 B spa- 
tial correlation. By integrating the expected spectra over 
4 < T vis < 20 MeV, e B ( 8 Li) and £b( 8 B) are estimated to be 
(77.6 ± 0.9)% and (88.4 ± 0.7)%, respectively. For show- 
ering muons, es is 100% because no correlation require- 
ment is imposed. For nonshowering muons, £5 is estimated 
from the 12 B and 12 N analysis (Sec. IV At , and the system- 
atic error arising from variations between isotopes is esti- 
mated with a FLUKA simulation. Figure [10] shows the im- 
pact parameter (AL) distribution for the 12 B and 12 N can- 
didate events for nonshowering muons (AC < 10 6 p.e.). We 
find that (91.6 ± 4.3)% of the candidates are within 3 m of 
the muon track. This fraction is the value of es( 12 B) for 
AL < 3 m. 

To obtain es( 8 Li) and es( 8 B), an additional correction for 
the difference between the muon- 8 Li or muon- 8 B and the 
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FIG. 1 1 : Allowed regions for the production rates of 8 Li and 8 B from 
the combined fits of the energy spectra and the AT distributions at 
la, 2a, and 3a C.L. for (a) showering muons (AC > 10 6 p.e.) and 
(b) nonshowering muons (AC < 10 B p.e.) with AL < 3 m track cut. 
The black points indicate the best-fit parameters. 
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muon- 12 B spatial correlations is applied. This correction is 
derived from the FLUKA simulation described in Sec. IVIBI 
The simulation does not include the uncertainties in the muon 
track and the isotope decay vertex reconstruction; it is only 
used to study the isotope dependance of eg- The range of 
values of e$ from FLUKA for different spallation isotopes is 
used to estimate the systematic error that should be added to 
es( 12 B) in order to obtain a common es for all spallation iso- 
topes. The resulting value, es — (91.6 ± 8.4)%, is used for 
estimating the 8 Li and 8 B (and, later the 9 C, 8 He, and 9 Li) 
yields. 

Combining the above analyses, we obtain V( 8 Li) = 
(12.2 ± 2.6) x 1(T 7 and Y( 8 B) = (8.4 ± 2.4) x 1CT 7 
/x _1 g _1 cm 2 . The isotope production rates are i?( 8 Li) = 
15.6±3.2andi?( 8 B) = 10.7±2.9kton~ 1 day- 1 . The contour 
plots in Fig.fTTIshow the correlation between 8 Li and 8 B. Due 
to their similar lifetimes, 8 Li and 8 B are identified primarily 
by their energy spectra. 

C. 8 He and 9 Li 

8 He (r= 171.7 ms, Q = 10.7 MeV) Etl and 9 Li 
(r = 257.2 ms, Q = 13.6 MeV) El /3"-decay candidate 
events are selected according to the cuts 1 < fJ V i S < 13 MeV 
and AT < 10 s, and by the detection of a neutron follow- 
ing the /3~-decay event. 8 He decays to neutron-unstable 
excited states of 8 Li with a (16 ± 1)% branching ratio 14111 . 
and 9 Li decays to neutron-unstable excited states of 9 Be 
with a (50.8 ± 0.9)% branching ratio |E1. The neutron is 
identified by the 2.225-MeV 7 ray from radiative capture 
on X U (1.8 < E vis < 2.6 MeV). The 7 ray is required to 
be within 200 cm and 1.0 ms of the 8 He or 9 Li /3~ -decay 
candidate. Finally, the 8 He- 9 Li analysis is performed using a 
5.5-m-radius spherical fiducial volume to reduce the number 
of accidental coincidences between the /3~-decay candidate 
and external 7 ray backgrounds near the balloon. 

The inset in Fig. Q~2] shows the AT distribution for the 
events that satisfy the criteria outlined above. Figure [T2l also 
shows the residual E v i s distribution corresponding to the sub- 



FIG. 12: Background-subtracted T v i s spectrum above 1 MeV, where 
signal and background events are taken from 0.002 < AT < 1 s 
and 5.002 < AT < 6 s, respectively. Detection of a neutron cap- 
ture following a /3-decay is required to select events from the 
/3 _ + n decay mode. The branching ratios for 9 Li and 8 He are 
(50.8 ± 0.9)% and (16 ± 1)%, respectively. The production rates 
of 8 He (t = 171.7 ms, Q = 10.7 MeV) and 9 Li (r = 257.2 ms, 
Q = 13.6 MeV) are estimated to be 1.0 ± 0.5 and 2.8 ± 0.2 
kton _1 day _1 , respectively, from simultaneously fitting the E v i B 
spectrum and the AT distribution shown in the inset. The fit to the 
AT distribution has a xV d -°- f ■= 95/97. 




FIG. 13: Allowed regions for the production rates of 9 Li and 8 He 
from the combined fits of the energy spectra and the AT distributions 
at la, 2a, and 3a C.L. for (a) showering muons (AC > 10 6 p.e.) and 
(b) nonshowering muons (AC < 10 6 p.e.) with AL < 3 m track cut. 
The black points indicate the best-fit parameters. 



traction of a background spectrum in the 5.002 < AT < 6 s 
window from a signal in the 0.002 < AT < 1 s window. The 
expected E v i s distributions for 8 He and 9 Li are calculated by 
incorporating the KamLAND response and adjusting for the 
energy deposited by the thermalizing neutron from 8 He or 9 Li 
decay. iV( 8 He) and iV( 9 Li) are determined from a simultane- 
ous binned maximum likelihood fit to the AT distribution and 
a chi-square fit to the E vis distribution. For the fit to the E vis 
distribution, the uncertainty in the energy scale parameters are 
treated in the same manner as the 8 Li- 8 B analysis described in 
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FIG. 14: Background-subtracted T v i s spectrum above 
12 MeV, where signal and background events are taken from 
0.2 < AT < 0.6 s and 10.2 < AT < 10.6 s, respectively. The 
production rate of 9 C (r = 182.5 ms, Q = 16.5 MeV) is estimated 
to be 3.8 ± 1.5 kton _1 day _1 from a simultaneous fit to the E v - m 
spectrum and the AT distribution shown in the inset. The fit to the 
AT distribution has a X 2 /d.o.f.= 97/96. 



FIG. 15: Background-subtracted T v j s spectrum above 1.4 MeV, 
where signal and background are selected from 5 < AT < 90 min 
and 185 < AT < 270 min, respectively. The production rate of 
U C (r = 29.4 min, Q = 1.98 MeV) is estimated to be 1106 ± 178 
kton _1 day _1 by fitting Eq. (6]l to the AT distribution shown in the 
inset. The fit to the AT distribution has a x 2 /d.o.f.= 76/57. 



Sec. ED 

The procedure for calculating the 8 He and 9 Li selec- 
tion efficiency is the same as for the 8 Li and 8 B effi- 
ciency analysis (Sec. IV BK except for the correction for 
the neutron detection requirement, which is calculated with 
the GEANT4-based Monte Carlo simulation described in 
Sec. lVTAl The resultant efficiencies e( 8 He) = (14.9 ± 1.0)% 
and e( 9 Li) = (46.1 ± 1.1)% include the appropriate f3-n 
branching fractions. Since a reduced volume is used 
in this analysis, the 1.6% fiducial volume uncertainty 
from Ref. 04411 is included in the above efficiencies. 
The resulting yields are y( 8 He) = (0.7 ± 0.4) x 10- y and 
F( 9 Li) = (2.2 ±0.2) x 10- 7 ^g^cm 2 . The produc- 
tion rates are i?( 8 He) = 1.0 ± 0.5 and i?( 9 Li) = 2.8 ± 0.2 
kton _1 day _1 . The contour plots in Fig. [13] show the corre- 
lation between 9 Li and 8 He. 



D. 9 C 

The inset in Fig. Q~4] shows the AT distribution for all 
events with visible energy 12 < E v i s < 20 MeV. The anal- 
ysis region (0.2 < AT < 0.6 s) contains events from 9 C 
(t = 182.5 ms, Q = 16.5 MeV) EH /?+ decay. iV( 9 C) is de- 
termined from a simultaneous binned maximum likelihood fit 
to the AT distribution and a chi-square fit to the E V [ S distri- 
bution. The uncertainty in the energy scale parameters are 
treated in the same manner as the 8 Li- 8 B analysis described 
in Sec. IV Bl In this fit, 8 Li, 8 B, and 9 Li are treated as pos- 
sible contaminants, the amounts are constrained to the val- 
ues obtained in the previously described analyses. This con- 
straint includes the correlation between 8 Li and 8 B shown 



in Fig. [TT| By integrating the theoretical 9 C E v - ls spec- 
trum, we obtain the efficiency for the 12 < E v i s < 20 MeV 
cut of e( 9 C) = (7.2 ± 1.0)%. Combining this with the above 
results gives y( 9 C) = (3.0 ± 1.2) x lO" 7 ^g^cm 2 and 
i?( 9 C) = 3.8 ± 1.5 kton-May- 1 . 

E. U C 

The production of U C ((3 + decay, r = 29.4 min, 
Q= 1.98 MeV) ll40Tl through muon-initiated spallation is 
usually accompanied by a neutron, allowing identification by 
the triple coincidence of the primary muon, the spallation neu- 
tron, and the subsequent j3 + 1 1 3l 14511 . The n C (3 + decays 
are selected in the range 1.4 < E V i S < 2.0 MeV and are pre- 
ceded by a detected muon that is accompanied by at least one 
neutron capture, identified by the 2.225-MeV 7 ray from cap- 
ture on 1 H. The 7 ray is required to be in the time window 
10 < AT < 2500^s relative to the muon. To reduce the back- 
ground, a 7-m-diameter fiducial volume is used. To avoid in- 
efficiencies from run boundaries and the long lifetime of n C, 
the first 5 h of the typically 24-h-long run are not used in the 
selection of the n C candidates. The number of muon- 11 C co- 
incidences is extracted from the AT distribution for all events 
that meet the criteria, shown in the inset of Fig. [TBI 

The efficiency determination takes into account the visible 
energy range for n C (3 + decay (22.7 ± 3.6)% and the previ- 
ously discussed neutron detection efficiency (Sec. |IV) . The 
efficiency also takes into account a correction for 11 C pro- 
duction modes, designated invisible modes, which do not pro- 
duce neutrons lfl3l l45tl . To measure this correction, muon- 
n C event pairs were selected with and without the neutron 
requirement for a subset of the data where the n C candidate 
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FIG. 16: Background-subtracted _B v i s spectrum above 2 MeV, 
where signal and background are selected using 10 < AT < 90 s 
and 190 < AT < 270 s, respectively. The production rate of 10 C 
(r = 27.8 s, Q = 3.65 MeV) is 21.1 ± 1.8 kton^day -1 , as deter- 
mined by fitting Eq. l|6} to the AT distribution shown in the inset. 
The fit to the AT distribution has a x 2 /d.o.f - 80/56. 
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FIG. 17: Background-subtracted T v i s spectrum above 5.5 MeV for 
showering muons (A£ > 10 6 p.e.), where signal and background 
are selected by 8 < AT < 60 s and 408 < AT < 460 s, respec- 
tively. The production rate of "Be (r = 19.9 s, Q = 11.5 MeV) 
for showering muons is estimated to be 1.0 ±0.2 kton _1 day _1 
by fitting Eq. {6} to the AT distribution shown in the inset 
(X 2 /cl.o.f.= 37/41). 



is required to be within 50 cm of the muon track; restricting 
the study to a subset of the data mitigated the reduced signal- 
to-background ratio associated with relaxing the neutron re- 
quirement. The number of muon- 11 C coincidences in each 
case was extracted from a fit of Eq. © to the corresponding 
AT distribution. The visible mode efficiency, e v j s , taken as 
the ratio of the number of muon- 11 C pairs with one or more 
neutrons to the number of muon- 11 C pairs without the neu- 
tron requirement, is (88.4 ± 2.4)%. Applying the correction 
for post-muon electronics effects, the visible mode fraction 
is (96.3 ± 2.0)%, consistent with Ref. [45], which obtains 
£vis = 95.6% for 285-GeV muons. 

Due to the relatively long n C lifetime, we also consid- 
ered the effect of diffusion. An analysis of 222 Rn that was 
accidentally introduced into the center of KamLAND dur- 
ing the deployment of a calibration device shows that the 
diffusion speed is approximately 1 mm/h. From this study, 
the effect of n C diffusion on the efficiency is estimated 
to be less than 0.5%. Combining this with the above re- 
sults gives Y( n C) = (866 ± 153) x 10~ 7 /Lt^g-W 2 and 
ft( n C) = 1106 ± 178 kton^day" 1 . 



F. 10 C 

As with n C, the production of 10 C(/3+ decay, r = 27.8 s, 
Q = 3.65 MeV) |4l|] through muon-initiated spallation is 
usually accompanied by a neutron, so the selection criterion 
requiring a triple coincidence of the primary muon, the neu- 
tron, and the 10 C candidate is used. The neutron is identi- 
fied by the 2.225-MeV n+ 1 H capture 7 ray. The number of 
10 C candidates is determined from fitting Eq. (O to the AT 



distribution for all events identified as 10 C, shown in the in- 
set in Fig. [16] 11 Be is a potential background for this 10 C 
analysis, but the correction to the 10 C yield is estimated to 
be less than 1% because of the low 11 Be production rate and 
the neutron coincidence requirement. The efficiency for the 
visible energy cut 2.0 < E vis < 4.0 MeV is (73.5 ± 3.2)%. 
The visible mode efficiency is e v ; s = (90.7 ± 5.5)% after cor- 
recting for the electronics effects following muons. The fi- 
nal efficiency is (89.6 ± 5.5)%. The resulting isotope yield is 
F( 10 C) = (16.5 ± 1.9) x 10~ 7 ^g^cm 2 and the produc- 
tion rate is i?( 10 C) = 21.1 ± l.Skton-May" 1 . 



G. n Be 

The u Be /?--decay (t = 19.9 s, Q = 11.5 MeV) Q 
events are selected according to 5.5 < E V i S < 16.0 MeV. The 
E v i s cut efficiency is 63.4%. The inset in Fig. [P71 shows the 
AT distributions for the events only after showering muons. 
For nonshowering muons, a tighter muon track cut AL < 1 m 
is applied in order to reduce the background rate. The track cut 
efficiency is estimated from the 12 B candidates using an anal- 
ysis similar to that in Sec. IVBI The resulting isotope yield is 
F( n Be) = (1.1 ± 0.2) x 10~ 7 ^g^cm 2 and production 
rate is i?( n Be) = 1.4 ± 0.3 kton-^ay- 1 . 



VI. MONTE CARLO SIMULATIONS 

The GEANT4 and FLUKA simulations are used to reproduce 
the measurements from KamLAND. While GEANT4 is used 
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FIG. 18: E v i s distribution of the prompt events from candidates iden- 
tified as neutrons produced by muon-induced spallation in the mate- 
rial outside of the KamLAND ID. 



only to simulate neutron production, both neutron and light 
isotope production are tested with FLUKA. 



FIG. 19: Radial distribution (from the center of KamLAND) 
of the prompt event from candidates identified as neutrons pro- 
duced by muon-induced spallation in the material outside of the 
KamLAND ID. Assuming an exponential distribution, a fit to the 
measured data (red dashed lines) yields an attenuation length of 
70 ± 2 g/cm 2 . A similar fit to the Monte Carlo (black dotted lines) 
yields 69 ±2 g/cm 2 . 



A. GEANT4 

GEANT4 is a widely used toolkit for performing particle 
tracking simulations on an event-by-event basis. A descrip- 
tion of the available physics processes included is given in 
Refs. H2IEII]. Here we compare the GEANT4 (version 9.1) 
prediction for neutron yield by spallation with the results ob- 
tained in Sec. [IV] We use the physics list QGS_BIC, de- 
veloped by the GEANT4 group to support the binary cas- 
cade (BIC) model at lower energies (below 10 GeV for p 
and n, and below 1.2 GeV for tt). This treatment is also 
appropriate for the simulation of interactions of nucleons 
and ions. At higher energies, a quark-gluon string (QGS) 
model is applied for the hadronic interactions. Neutron elas- 
tic and inelastic interactions below 20 MeV are described 
by a high-precision data-driven model (NeutronHP). The 
G4EMEXTRAPHYSICS physics list is also used to model the 
photonuclear and muon-nuclear interaction processes, which 
dominate the neutron production by muons in the simulation. 

To estimate the neutron production yield as a function of 
muon energy in GEANT4, monoenergetic muons of several en- 
ergies are injected at the center of a generic hydrocarbon block 
of thickness 40 m. The region more than 10 m away from the 
edges of the block is analyzed to avoid boundary effects. As 
shown later in Fig. [20] the neutron production yields predicted 
by GEANT4 are systematically lower than experiment, with 
the exception of point (F) from the LVD lll2ll . These results 
are consistent with previous work [46]. 

A Monte Carlo simulation based upon GEANT4 and MU- 
SIC (described in Sec. ITTTb is used to study neutrons produced 
by muon-induced spallation in the material outside of the 
KamLAND ID. Some of these neutrons have sufficient energy 



to enter the ID where they thermalize and capture. They can 
be identified by the coincidence of a prompt signal (for exam- 
ple from n + p elastic scattering) and a delayed signal from 
the capture 7 ray. This is the same inverse /3-decay reaction 
signature used for T7 e detection, V e + p — > e + + n, where the 
e + is the prompt signal, and the 7 ray from neutron capture 
is the delayed signal and therefore a potential background. 
These neutrons are also a background for dark-matter experi- 
ments that employ nuclear recoils as a detection method. 

The primary purpose of this Monte Carlo simulation is 
to estimate the rate of untagged fast neutrons, i.e. neutrons 
produced by muons where the muon is undetected by the 
KamLAND ID. A few of these muons are detected by the OD, 
either from the Cherenkov radiation produced by the muon 
itself (tracked muons), or by accompanying electromagnetic 
and hadronic showers that enter the OD (untracked muons). A 
prompt signal in coincidence with these tracked and untracked 
muons followed by a delayed capture 7 signal identifies can- 
didates for untagged fast neutrons. 

Fast neutrons generated by untracked muons are produced 
primarily in the rock surrounding the OD, whereas fast neu- 
trons from tracked muons are primarily produced in the water 
of the OD. It is shown in the Monte Carlo that the tracked and 
untracked muons give distinguishable OD visible energy dis- 
tributions. A comparison between Monte Carlo and measure- 
ment of the distribution of the number of OD PMTs with sig- 
nals above threshold for tracked and untracked muons reveals 
a deficiency of fast neutrons from untracked muons, consistent 
with the underproduction of neutrons by GEANT4 in concrete 
reported in Ref. [47]. 

Figure [18] shows the E v - ls distribution of the prompt events 
from the Monte Carlo simulation of tracked and untracked 
muons compared with the measured data. The measured data 
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FIG. 20: Neutron production yield in liquid scintillator as a function 
of muon energy. The red star point shows the KamLAND result, 
Y n = (2.8 ± 0.3) x 10- 4 ^- 1 g _1 cm 2 , for 260 ± 8 GeV. Other 
points show the results from experiments at (A) 20 mwe 1481 . (B) 
25 mwe (U, (C) 32 mwe (U, (D) 316 mwejll, (E) 570 mwe 
(U, (F) 3000 mwe d, and (G) 5200 mwe fullThe blue square 
and green triangle show the GEANT4 and FLUKA Monte Carlo pre- 
dictions, respectively, from monochromatic muon beams. 



and the Monte Carlo simulation correspond to an equal live 
time exposure of 1368 days. Figure [T9lshows the radial distri- 
bution of the prompt events relative to the KamLAND center. 
Both the measured data and the GEANT4 simulation exhibit an 
exponential attenuation of the neutrons as they penetrate far- 
ther into the detector, with the simulation yielding an attenua- 
tion length of 69 ± 2 g/cm 2 , consistent with the measurement 
70 ± 2 g/cm 2 . 



B. FLUKA 

FLUKA is a mature code that models nuclear and particle 
physics processes from thermal neutrons to heavy-ion colli- 
sions lfl5l [HI] . It has been used previously to model muon- 
initiated spallation in liquid scintillator HH EHH HH . We 
use FLUKA version 2006.3b to model neutron and light isotope 
production from muon-initiated spallation in KamLAND. A 
40-m-radius by 40-m-high cylinder of KamLAND liquid scin- 
tillator is used in the simulation; the concentric inner cylinder 
of 20-m radius and 20-m length is used for analysis. 

To estimate neutron production yield as a function of muon 
energy in FLUKA, monoenergetic beams of yT r anging from 
10 to 350 GeV were simulated as in Refs. JU H |H HI . 
Care is taken not to double-count the neutrons involved in 
reactions like (n, In). The results of this simulation are in- 



cluded in Fig. [20] The neutron production yield of this FLUKA 
simulation is 10% lower than previous work |[3Tll46ll52T,l53ll . 
but the power-law dependance on muon energy (E™, where 
a = 0.77) is consistent. Different scintillator compositions 
were studied, but they could not explain the deficit. This 
deficit is insignificant compared to the discrepancy between 
these simulations and the data. 

The production of light isotopes was studied in the same 
simulation. The results, including the primary production pro- 
cess and power-law exponent, are summarized in Table ITVl 
For some isotopes the primary production process is much 
larger than any secondary processes, as in the case of 12 B. For 
other isotopes the primary production process is only slightly 
larger than the secondary processes, as is the case for 9 Li. The 
isotopes produced primarily by 7 interactions, n C and 10 C, 
show the weakest dependence on muon energy. In compari- 
son, 12 N and 13 N, where the primary production mechanism 
is by p interactions, show the strongest dependence on muon 
energy. 

The use of a monoenergetic beam overestimates the 
production of neutrons and light isotopes. Simulations us- 
ing a monoenergetic /i + beam and a beam with the energy 
spectrum from Ref. 113311 were also run. The simulations show 
that the production yield for n + relative to /i~ is on average 
0.96±0.01 for the light isotopes. This reduction is expected, 
since /i~ may capture, creating spallation products, while fi + 
may not. This ratio, combined with the fi + to /i~ ratio at 
KamLAND, leads to a correction to the flux of 0.98 ±0.06 
for light isotopes and 0.981 ±0.005 for neutrons. The reduced 
production yield due to averaging over the muon spectrum is 
on average 0.92±0.02 for the light isotopes, which is slightly 
higher than the correction factor suggested by Ref. jlOll . The 
results of the FLUKA simulations for KamLAND presented in 
Table [Vlinclude these two corrections. 



VII. DISCUSSION 

The isotope production yields from muon-initiated spalla- 
tion in a liquid-scintillator target were investigated at CERN 
by earlier experiment IToll using the SPS muon beam with 
muon energies of 100 and 190 GeV. Based on those cross- 
section measurements and the predicted muon energy spec- 
trum at the KamLAND site, we calculated the isotope pro- 
duction rates by extrapolation, assuming a power-law of the 
muon energy E^. The mean muon energy at KamLAND is 
260 ± 8 GeV, so KamLAND provides data to test the extrap- 
olation method at this energy. 

The production yields for the isotopes from muon-initiated 
spallation in KamLAND are provided in Table [V] On av- 
erage, the yields from the showering muons (~15% of all 
muons), whose excess light yield parameter [AC, Eq. (f2]l] is 
greater than 10 6 p.e. (-3 GeV), constitute (70 ± 2)% of the 
yield from all muons. The production yield for n C is the 
largest, and its measured yield is larger than the FLUKA cal- 
culation by a factor of ~2. The Borexino Collaboration also 
reported a similar discrepancy [2, 3], which is consistent with 
what is observed in KamLAND. Some measured production 
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TABLE IV: Simulation of neutron and light isotope production in the KamLAND LS by muon-initiated spallation with FLUKA: the isotope 
production yield by monoenergetic (260 GeV) the ratio of the production yields by monoenergetic ji + compared to fi~; the ratio of the 
production yields by a fj,~ spectrum that matches Fig. 6 (KamLAND curve) in Ref. [33] compared with monoenergetic \x~\ the power-law 
exponent for the production yield Y(E I1 ) oc E" from a fit to the yields from monoenergetic \i~ with 10 < E^ < 350 GeV; and the primary 
process for producing the isotope. The uncertainties are statistical. 



Simulated production yield Ratio of Simulated Production Yields for 





(xHTV^e^cm 2 ) 




Spectrum/monoenergetic 


Power-law exp. 


Primary process 


n 


2344 ± 4 


0.969 ± 0.002 


0.912 ± 0.003 


0.779 ± 0.001 


7T _ ± 1 H, 12 C 


u c 


460.8 ± 1.7 


0.971 ± 0.005 


0.913 ± 0.006 


0.703 ± 0.002 


12 C(7,n) 


7 Be 


116.8 ± 0.9 


0.986 ± 0.011 


0.945 ± 0.011 


0.684 ± 0.004 


12 C(7, no) 


10 Be 


44.63 ± 0.53 


0.960 ± 0.018 


0.891 ± 0.019 


0.825 ± 0.007 


12 C(n, 3 He) 


12g 


30.85 ± 0.44 


0.970 ± 0.021 


0.936 ± 0.022 


0.828 ± 0.009 


12 C(n,p) 


8 Li 


23.42 ± 0.39 


0.927 ± 0.026 


0.936 ± 0.025 


0.821 ± 0.010 


12 C(n,pa) 


io c 


21.13 ± 0.37 


0.982 ± 0.025 


0.915 ± 0.027 


0.810 ± 0.010 


V ' "f) 


6 He 


13.40 ± 0.29 


0.916 ±0.035 


0.918 ±0.035 


0.818 ±0.013 


12 C(n,2p 3 He) 


8 B 


6.40 ± 0.20 


0.996 ± 0.045 


0.915 ±0.050 


0.804 ±0.019 


12 C(tt + , 2 H 2 H) 


°Li 


3.51 ±0.15 


0.856 ± 0.074 


0.842 ± 0.078 


0.801 ±0.026 


12 C(7r-, 3 He) 


9 C 


1.49 ± 0.10 


0.850 ±0.114 


0.949 ±0.102 


0.772 ± 0.039 


12 C>+, 3 H) 


12 N 


0.86 ± 0.07 


0.963 ±0.128 


1.006 ±0.120 


0.921 ±0.045 


12 C(p,n) 


n Be 


0.94 ± 0.08 


0.842 ±0.145 


0.804 ±0.161 


0.753 ±0.051 


12 C(n,2p) 


8 He 


0.35 ± 0.05 


0.964 ± 0.200 


0.576 ± 0.372 


0.926 ± 0.078 


12 C(7r-,n3p) 


13g 


0.31 ±0.04 


1.020 ± 0.197 


1.062 ±0.176 


0.742 ± 0.075 


13 C(n,p) 


!5 


0.05 ± 0.02 


1.250 ± 0.379 


1.635 ±0.234 


0.793 ± 0.244 


16 0( 7 ,n) 


13 N 


0.06 ± 0.02 


1.500 ± 0.272 


1.190 ±0.401 


1.120 ±0.220 


13 C(p,n) 



TABLE V: Summary of the neutron and isotope production yields from muon-initiated spallation in KamLAND. The results of the FLUKA 
calculation include corrections for the muon spectrum and the /i + composition of the cosmic-ray muon flux. 





Lifetime in 


Radiation energy 




Yield (xlO" 7 ^"^ 




Fraction from showering /i 




KamLAND LS 


(MeV) 


Ref. QO] 


FLUKA calc. 


This measurement 


This measurement 




207.5 ^s 


2.225 (capt. 7) 




2097 ± 13 


2787 ± 311 


(64 ± 5)% 


12g 


29.1 ms 


13.4 08") 




27.8 ± 1.9 


42.9 ±3.3 


(68 ± 2)% 


12 N 


15.9 ms 


17.3 (/3+) 




0.77 ±0.08 


1.8 ± 0.4 


(77 ± 14)% 


8 Li 


1.21 s 


16.0 (j3-a) 


1.9 ±0.8 


21.1 ± 1.4 


12.2 ±2.6 


(65 ± 17)% 


8 B 


1.11 s 


18.0 (fi+a) 


3.3 ±1.0 


5.77 ±0.42 


8.4 ±2.4 


(78 ± 23)% 


"C 


182.5 ms 


16.5 08+) 


2.3 ±0.9 


1.35 ±0.12 


3.0 ± 1.2 


(91 ± 32)% 


8 He 


171.7 ms 


10.7 (/3~7n) 




0.32 ±0.05 


0.7 ±0.4 


(76 ± 45)% 


(, Li 


257.2 ms 


13.6 (/3~"/n) 


| 1.0 ±0.3 


3.16 ±0.25 


2.2 ± 0.2 


(77 ± 6)% 


U C 


29.4 min 


1.98 03+) 


421 ± 68 


416 ± 27 


866 ± 153 


(62 ± 10)% 


io c 


27.8 s 


3.65 03+7) 


54 ± 12 


19.1 ± 1.3 


16.5 ±1.9 


(76 ± 6)% 


n Be 


19.9 s 


11.5 on 


< 1.1 


0.84 ±0.09 


1.1 ±0.2 


(74 ± 12)% 


6 He 


1.16 s 


3.51 OS") 


7.5 ±1.5 


12.08 ±0.83 






7 Be 


76.9 day 


0.478 (EC 7) 


107 ± 21 


105.3 ±6.9 







yields, such as 8 Li and °C, deviate significantly from esti- 
mates based on the muon beam experiment, indicating that 
perhaps estimation by extrapolation is not sufficient. All iso- 
tope yields are consistent within an order of magnitude. 



VIII. SUMMARY 

We have analyzed KamLAND data to measure pro- 
duction yields of radioactive isotopes and neutrons 
through muon-initiated spallation in liquid-scintillator. 
The neutron production yield is evaluated to be 
Y n = (2.8 ± 0.3) x KTV^g^cm 2 , which is higher 
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than the expectation from Monte Carlo simulations based 
on GEANT4 and FLUKA. Some isotope production yields 
are found to be inconsistent with extrapolations — based on 
a power-law dependance with respect to muon energy — of 
results from muon beam experiments. 

Acknowledgments 

The KamLAND experiment is supported by the COE pro- 
gram under Grant 09CE2003 of the Japanese Ministry of Ed- 



ucation, Culture, Sports, Science and Technology; the World 
Premier International Research Center Initiative (WPI Ini- 
tiative), MEXT, Japan; and under the U.S. Department of 
Energy (DOE) Grants DEFG03-00ER41 138 and DE-AC02- 
05CH1 1231, as well as other DOE Grants to individual insti- 
tutions. We are grateful to the Kamioka Mining and Smelting 
Company for supporting the activities in the mine. We would 
also like to thank V. A. Kudryavtsev and M. Spurio for helpful 
comments. 



[1] G. Bellini, J. Benziger, S. Bonetti, M. B. Avanzini, B. Cac- 
cianiga, L. Cadonati, F. Calaprice, C. Carraro, A. Chavar- 
ria, F. Dalnoki-Veress, et al. (Borexino Collaboration), astro- 
ph/0808.2868vl. 

[2] C. Arpesella, G. Bellini, J. Benziger, S. Bonetti, B. Caccian- 
iga, F. Calaprice, F. Dalnoki-Veress, D. D' Angelo, H. de Kerret, 
A. Derbin, et al. (Borexino Collaboration), Phys. Lett. B 658, 
101 (2008). 

[3] C. Arpesella, H. O. Back, M. Balata, G. Bellini, J. Benziger, 
S. Bonetti, A. Brigatti, B. Caccianiga, L. Cadonati, F. Calaprice, 
et al. (Borexino Collaboration), Phys. Rev. Lett. 101, 091302 
(2008). 

[4] Y. Hirano, T. Kishimoto, I. Ogawa, R. Hazama, S. Umehara, 
K. Matsuoka, G. Ito, and Y. Tsubota, J. Phys. Conf. Ser. 120, 
052053 (2008). 

[5] T. Kishimoto, talk at the International Workshop on "Double 

Beta Decay and Neutrinos", June 11-13, 2007, Osaka, available 

online at http://dbd07.phys.sci.osaka-u.ac.jp/ (2007). 
[6] K. Zuber, Am. Inst. Phys. Conf. Proc. 942, 101 (2007). 
[7] M. C. Chen, Earth, Moon, and Planets 99, 221 (2006). 
[8] C. Grieb, J. M. Link, M. L. Pitt, R. S. Raghavan, D. Rountree, 

and R. B. Vogelaar, hep-ex/0705.2769vl. 
[9] T. M. Undagoitia, F. von Feilitzsch, M. Goger-Neff, L. Ober- 

auer, W. Potzel, A. Ulrich, J. Winter, and M. Wurm, J. Phys. 

Conf. Ser. 120, 052018 (2008). 
[10] T. Hagner, R. von Hentig, B. Heisinger, L. Oberauer, 

S. Schonert, F. von Feilitzsch, and E. Nolte, Astropart. Phys. 

14,33 (2000). 

[11] M. Aglietta et al., Nuovo Cimento Soc. Ital. Fis. C 12, 467 
(1989). 

[12] M. Aglietta, E. Alyea, P. Antonioli, G. Badino, G. Bari, 
M. Basile, V. Berezinsky, F. Bersani, M. Bertaina, R. Bertoni, 
et al., Phys. Atomic Nuclei 66, 123 (2003). 

[13] H. Back, M. Balata, G. Bellini, J. Benziger, S. Bonetti, B. Cac- 
cianiga, F. Calaprice, D. D'Angelo, A. de Bellefon, H. de Ker- 
ret, et al. (Borexino Collaboration), Phys. Rev. C 74, 045805 
(2006). 

[14] P. Antonioli, C. Ghetti, E. V. Korolkova, V. A. Kudryavtsev, and 

G. Sartorelli, Astropart. Phys. 7, 357 (1997). 
[15] A. Fasso, A. Ferrari, S. Roesler, P. Sala, G. Battistoni, F. Cerutti, 

E. Gadioli, M. Garzelli, F. Ballarini, A. Ottolenghi, et al., hep- 

ph/0306267. 

[16] A. Ferrari, P. R. Sala, A. Fasso, and J. Ranft, FLUKA: A 
multi-particle transport code (program version 2005) (CERN, 
Geneva, 2005). 

[17] S. Agostinelli, J. Allison, K. Amako, J. Apostolakis, H. Araujo, 
P. Arce, M. Asai, D. Axen, S. Banerjee, G. Barrand, et al., Nucl. 



Instr. Meth. A 506, 250 (2003). 
[18] J. Allison, K. Amako, J. Apostolakis, H. Araujo, 
P. Arce Dubois, M. Asai, G. Barrand, R. Capra, S. Chau- 
vie, R. Chytracek, et al., IEEE Trans. Nucl. Sci. 53, 270 
(2006). 

[19] R. Battino, T. R. Rettich, and T. Tominaga, J. Phys. Chem. Ref. 

Data 13, 563 (1984). 
[20] P. J. Hesse, R. Battino, P. Scharlin, and E. Wilhelm, J. Chem. 

Eng. Data 41, 195 (1996). 
[21] H. Kume, S. Sawaki, M. Ito, K. Arisaka, T. Kajita, 

A. Nishimura, and A. Suzuki, Nucl. Instr. Meth. Phys. Res. 205, 

443 (1983). 

[22] S. Kleinfelder, IEEE Trans. Nucl. Sci. 50, 955 (Aug. 2003). 

[23] D. W. McKee, J. K. Busenitz, and I. Ostrovskiy, Nucl. Instr. 
Meth. A 587, 272 (2008). 

[24] B. E. Berger, J. Busenitz, T. Classen, M. P. Decowski, D. A. 
Dwyer, G. Elor, A. Frank, S. J. Freedman, B. K. Fujikawa, 
M. Galloway, et al. (KamLAND Collaboration), JINST 4, 
P04017 (2009). 

[25] J. B. Birks, Proc. Phys. Soc. A64, 874 (1951). 

[26] J. B. Birks, The Theory and Practice of Scintillation Counting 
(Pergamon, London, 1964). 

[27] Digital Map 50 m Grid (Elevation), Geographical Survey Insti- 
tute of Japan (1997), unpublished. 

[28] Y. Becherini, A. Margiotta, M. Sioli, and M. Spurio, Astropart. 
Phys. 25, 1 (2006). 

[29] M. Honda, T. Kajita, K. Kasahara, S. Midorikawa, and 
T. Sanuki, Phys. Rev. D 75, 043006 (2007). 

[30] D. Casper, Nucl. Phys. B Proc. Suppl. 112, 161 (2002). 

[31] D.-M. Mei and A. Hime, Phys. Rev. D 73, 053004 (2006). 

[32] C. Galbiati and J. F. Beacom, Phys. Rev. C 72, 025807 (2005). 

[33] A. Tang, G. Horton-Smith, V. A. Kudryavtsev, and A. Tonazzo, 
Phys. Rev. D 74, 053007 (2006). 

[34] S. Eidelman, K. Hayes, K. Olive, M. Aguilar-Benitez, C. Am- 
sler, D. Asner, K. Babu, R. Barnett, J. Beringer, P. Burchat, 
et al., Phys. Lett. B 592, 1 (2004). 

[35] Tech. Rep., Kamioka Mining & Smelting Company (1977), in- 
ternal report. 

[36] D. E. Groom, N. V. Mokhov, and S. I. Striganov, Atomic Data 

and Nuclear Data Tables 78, 183 (2001). 
[37] P. H. Barrett, L. M. Bollinger, G. Cocconi, Y. Eisenberg, and 

K. Greisen, Rev. Mod. Phys. 24, 133 (1952). 
[38] S. F. Mughabghab, M. Divadeenam, and N. E. Holden, Neutron 

Cross Sections, Volume 1, Neutron Resonance Parameter and 

Thermal Cross Sections, Part A Z = 1 — 60 (Academic Press, 

New York, 1981). 
[39] S. Baker and R. D. Cousins, Nucl. Instr. Meth. Phys. Res. 221, 



16 



437 (1984). 

[40] F. Ajzenberg-Selove, Nucl. Phys. A 506, 1 (1990). 

[41] D. R. Tilley, J. H. Kelley, J. L. Godwin, D. J. Millener, J. E. 

Purcell, C. G. Sheu, and H. R. Weller, Nucl. Phys. A 745, 155 

(2004). 

[42] W. T. Winter, S. J. Freedman, K. E. Rehm, and J. P. Schiffer, 

Phys. Rev. C 73, 025503 (2006). 
[43] M. Bhattacharya, E. G. Adelberger, and H. E. Swanson, Phys. 

Rev. C 73, 055802 (2006). 
[44] S. Abe, T. Ebihara, S. Enomoto, K. Furuno, Y. Gando, 

K. Ichimura, H. Ikeda, K. Inoue, Y. Kibe, Y. Kishimoto, 

et al. (KamLAND Collaboration), Phys. Rev. Lett. 100, 221803 

(2008). 

[45] C. Galbiati, A. Pocar, D. Franco, A. Ianni, L. Cadonati, and 

S. Schonert, Phys. Rev. C 71, 055805 (2005). 
[46] H. Araujo, V. Kudryavtsev, N. Spooner, and T. Sumner, Nucl. 

Instr. Meth. A 545, 398 (2005). 
[47] M. G. Marino, J. A. Detwiler, R. Henning, R. A. Johnson, A. G. 



Schubert, and J. F. Wilkerson, Nucl. Instr. Meth. A 582, 611 
(2007). 

[48] R. Hertenberger, M. Chen, and B. L. Dougherty, Phys. Rev. C 
52,3449 (1995). 

[49] L. B. Bezrukov, V. I. Beresnev, G. T. Zatsepin, O. G. Ryazh- 

skaya, and L. N. Stepanets, Sov. J. Nucl. Phys. 17, 51 (1973). 
[50] F. Boehm, J. Busenitz, B. Cook, G. Gratta, H. Henrikson, J. Ko- 

rnis, D. Lawrence, K. B. Lee, K. McKinny, L. Miller, et al., 

Phys. Rev. D 62, 092005 (2000). 
[51] R. I. Enikeev, G. T. Zatsepin, E. V. Korol'kova, V. A. 

Kudryavstev, A. S. Malgin, O. G. Ryazhskaya, and F. F. 

Khal'chukov, Sov. J. Nucl. Phys. 46, 883 (1987). 
[52] Y.-F. Wang, V. Balic, G. Gratta, A. Fasso, S. Roesler, and 

A. Ferrari, Phys. Rev. D 64, 013012 (2001). 
[53] V. A. Kudryavtsev, N. J. C. Spooner, and J. E. McMillan, Nucl. 

Instr. Meth. A 505, 688 (2003). 



